clear all
set more off
*cd /Users/zimaoxiao/Dropbox/CC_Yield_Predict/replication_package/  // Navigate to replication folder on your own machine
cd /Users/zimaoxiao/Desktop/replication_package
global weather "timeInt0_3 timeInt3_6 timeInt6_9 timeInt9_12 timeInt12_15 timeInt15_18 timeInt18_21 timeInt21_24 timeInt24_27 timeInt27_30 timeInt30_33 timeInt33_36 timeInt36_39 prec prec2"

qui {
* ============================
* 	  Linear trends spec.
* ============================
use if inrange(year,1950,2020) using data/regression/reg_2050_prediction_3Cstep, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
	
* regression
reghdfe lncornyield $weather state_fips#c.trend [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_lin
parmest, label saving(`"`corn_st_lin'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)

* extract temperature coefficients		
use parm estimate using `corn_st_lin', clear
keep if _n<=13
forvalues i = 0(3)36 {
	local j = `i' + 3
	gen timeInt`i'_`j' = estimate if parm == "timeInt`i'_`j'"
	sum timeInt`i'_`j'
	local k = r(mean)
	replace timeInt`i'_`j' = `k'
}
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace

* extract county linear trends coefficients		
use parm estimate using `corn_st_lin', clear		
drop if _n<16
gen quad = ustrright(parm,1)
preserve
drop if quad != "d"
drop quad
destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
drop parm			
rename estimate stlin
tempfile cornlintrend
save `cornlintrend', replace
restore
						
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2050_prediction_3Cstep, clear
bys county_fips: keep if _n==1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
cross using `corntempcoef'

foreach k in ssp245 ssp585 {
	preserve
	* calculate climate change impacts by county
	bys county_fips: gen ccimpact_corn_`k' = (CC_timeInt0_3_`k'*timeInt0_3) + (CC_timeInt3_6_`k'*timeInt3_6) + (CC_timeInt6_9_`k'*timeInt6_9) + (CC_timeInt9_12_`k'*timeInt9_12) + (CC_timeInt12_15_`k'*timeInt12_15) + (CC_timeInt15_18_`k'*timeInt15_18) + (CC_timeInt18_21_`k'*timeInt18_21) + (CC_timeInt21_24_`k'*timeInt21_24) + (CC_timeInt24_27_`k'*timeInt24_27) + (CC_timeInt27_30_`k'*timeInt27_30) + (CC_timeInt30_33_`k'*timeInt30_33) + (CC_timeInt33_36_`k'*timeInt33_36) + (CC_timeInt36_39_`k'*timeInt36_39)
					
	* calculate overall impacts by county (climate change + tech effects)
	* 1950 as the 1st year, t=1; 2020 as the 71st year, t=71; 2050 as the 101th year, t=101
	bys county_fips: gen trends_corn_`k' = (CC_timeInt0_3_`k'*timeInt0_3) + (CC_timeInt3_6_`k'*timeInt3_6) + (CC_timeInt6_9_`k'*timeInt6_9) + (CC_timeInt9_12_`k'*timeInt9_12) + (CC_timeInt12_15_`k'*timeInt12_15) + (CC_timeInt15_18_`k'*timeInt15_18) + (CC_timeInt18_21_`k'*timeInt18_21) + (CC_timeInt21_24_`k'*timeInt21_24) + (CC_timeInt24_27_`k'*timeInt24_27) + (CC_timeInt27_30_`k'*timeInt27_30) + (CC_timeInt30_33_`k'*timeInt30_33) + (CC_timeInt33_36_`k'*timeInt33_36) + (CC_timeInt36_39_`k'*timeInt36_39) + (stlin*30)
	
	keep county_fips ccimpact_corn_`k' trends_corn_`k' 	
	tempfile corn_cccoef_`k'
	save `corn_cccoef_`k'', replace
	restore	
}
	
* merge in regression coefficients, calculate predicted 2020 yields
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef_ssp245', keep(3) noreport nogenerate
merge 1:1 county_fips using `corn_cccoef_ssp585', keep(3) noreport nogenerate
	
* ssp245
* predicted 2050 yields: climate change
bys county_fips: gen cc2050_lin_245 = cornyield_2020*(1+ccimpact_corn_ssp245)

* predicted 2050 yields: climate change + tech effects
bys county_fips: gen trends2050_lin_245 = cornyield_2020*(1+trends_corn_ssp245)

* ssp585
* predicted 2050 yields: climate change
bys county_fips: gen cc2050_lin_585 = cornyield_2020*(1+ccimpact_corn_ssp585)

* predicted 2050 yields: climate change + tech effects
bys county_fips: gen trends2050_lin_585 = cornyield_2020*(1+trends_corn_ssp585)
		
* Table S3: rows (1), (2), (4), (6), (8)
rename (cornyield_2020 ccimpact_corn_ssp245 trends_corn_ssp245 ccimpact_corn_ssp585 trends_corn_ssp585) (yield2020 ccimpact_corn_245_lin trends_corn_245_lin ccimpact_corn_585_lin trends_corn_585_lin)
tempfile baseline_2050_lin
save `baseline_2050_lin', replace
		
* ============================
* 	 Quadratic trends spec.
* ============================
use if inrange(year,1950,2020) using data/regression/reg_2050_prediction_3Cstep, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
				
* regression
reghdfe lncornyield $weather state_fips#c.trend state_fips#c.trend2 [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_quad
parmest, label saving(`"`corn_st_quad'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)
											
* extract temperature coefficients		
use parm estimate using `corn_st_quad', clear
keep if _n<=13
forvalues i = 0(3)36 {
	local j = `i' + 3
	gen timeInt`i'_`j' = estimate if parm == "timeInt`i'_`j'"
	sum timeInt`i'_`j'
	local k = r(mean)
	replace timeInt`i'_`j' = `k'
}
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace
						
* extract county linear trends coefficients		
use parm estimate using `corn_st_quad', clear
drop if _n<16
gen quad = ustrright(parm,1)
preserve
	* first-order
	drop if quad != "d"
	drop quad
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	drop parm			
	rename estimate stlin
	tempfile cornlintrend
	save `cornlintrend', replace
restore
	* second-order
	keep if quad == "2"
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	gen state_fips_1 = floor(state_fips/10)
	replace state_fips = state_fips_1
	drop state_fips_1
	drop parm quad
	rename estimate stquad
	tempfile cornquadtrend
	save `cornquadtrend', replace
	
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2050_prediction_3Cstep, clear
bys county_fips: keep if _n==1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
merge m:1 state_fips using `cornquadtrend', assert(3) keep(3) noreport nogenerate
cross using `corntempcoef'
	
foreach k in ssp245 ssp585 {
	preserve
	* calculate climate change impacts by county
	bys county_fips: gen ccimpact_corn_`k' = (CC_timeInt0_3_`k'*timeInt0_3) + (CC_timeInt3_6_`k'*timeInt3_6) + (CC_timeInt6_9_`k'*timeInt6_9) + (CC_timeInt9_12_`k'*timeInt9_12) + (CC_timeInt12_15_`k'*timeInt12_15) + (CC_timeInt15_18_`k'*timeInt15_18) + (CC_timeInt18_21_`k'*timeInt18_21) + (CC_timeInt21_24_`k'*timeInt21_24) + (CC_timeInt24_27_`k'*timeInt24_27) + (CC_timeInt27_30_`k'*timeInt27_30) + (CC_timeInt30_33_`k'*timeInt30_33) + (CC_timeInt33_36_`k'*timeInt33_36) + (CC_timeInt36_39_`k'*timeInt36_39)
		
	* calculate overall impacts by county (climate change + tech effects)
	* 1950 as the 1st year, t=1; 2020 as the 71st year, t=71; 2050 as the 101th year, t=101
	bys county_fips: gen trends_corn_`k' = (CC_timeInt0_3_`k'*timeInt0_3) + (CC_timeInt3_6_`k'*timeInt3_6) + (CC_timeInt6_9_`k'*timeInt6_9) + (CC_timeInt9_12_`k'*timeInt9_12) + (CC_timeInt12_15_`k'*timeInt12_15) + (CC_timeInt15_18_`k'*timeInt15_18) + (CC_timeInt18_21_`k'*timeInt18_21) + (CC_timeInt21_24_`k'*timeInt21_24) + (CC_timeInt24_27_`k'*timeInt24_27) + (CC_timeInt27_30_`k'*timeInt27_30) + (CC_timeInt30_33_`k'*timeInt30_33) + (CC_timeInt33_36_`k'*timeInt33_36) + (CC_timeInt36_39_`k'*timeInt36_39) + (stlin*30) + (stquad*5160)

	keep county_fips ccimpact_corn_`k' trends_corn_`k' 
	tempfile corn_cccoef_`k'
	save `corn_cccoef_`k'', replace
	restore	
}
	
* merge in regression coefficients, calculate predicted 2050 yields
* 1
* calculate predicted 2050 yield
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef_ssp245', keep(3) noreport nogenerate
merge 1:1 county_fips using `corn_cccoef_ssp585', keep(3) noreport nogenerate
	
* ssp245
* predicted 2050 yields: under climate change
bys county_fips: gen cc2050_quad_245 = cornyield_2020*(1+ccimpact_corn_ssp245)

* predicted 2050 yields: under climate change + tech effects
bys county_fips: gen trends2050_quad_245 = cornyield_2020*(1+trends_corn_ssp245)

* ssp585
* predicted 2050 yields: under climate change
bys county_fips: gen cc2050_quad_585 = cornyield_2020*(1+ccimpact_corn_ssp585)

* predicted 2050 yields: under climate change + tech effects
bys county_fips: gen trends2050_quad_585 = cornyield_2020*(1+trends_corn_ssp585)

* Table S3: rows (3), (5), (7), (9)	
rename (cornyield_2020 ccimpact_corn_ssp245 trends_corn_ssp245 ccimpact_corn_ssp585 trends_corn_ssp585) (yield2020 ccimpact_corn_245_quad trends_corn_245_quad ccimpact_corn_585_quad trends_corn_585_quad)
tempfile baseline_2050_quad
save `baseline_2050_quad', replace

* ============================
* 	         Output
* ============================
use `baseline_2050_lin', clear
merge 1:1 county_fips using `baseline_2050_quad', assert(3) keep(3) nogenerate
noi logout, save(output/TableS3) word dec(1) replace: tabstat yield2020 cc2050_lin_245 cc2050_quad_245 trends2050_lin_245 trends2050_quad_245 cc2050_lin_585 cc2050_quad_585 trends2050_lin_585 trends2050_quad_585, stat(p10 p25 p50 p75 p90) format(%10.4f) c(s)
save data/figure/baseline_2050, replace

}
*** EOF
